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Abstract 

This paper explores the relationship that exists between nonlinear normal modes (NNMs) 
defined as invariant manifolds in phase space and the spectral expansion of the Koopman 
operator. Specifically, we demonstrate that NNMs correspond to zero level sets of specific 
eigenfunctions of the Koopman operator. Thanks to this direct connection, a new, global 
parametrization of the invariant manifolds is established. Unlike the classical parametriza- 
tion using a pair of state-space variables, this parametrization remains valid whenever the 
invariant manifold undergoes folding, which extends the computation of NNMs to regimes 
of greater energy. The proposed ideas are illustrated using a two-degree-of-freedom system 
with cubic nonlinearity. 


1. Introduction 

Nonlinear normal modes (NNMs) of vibration can be regarded as the extension of the 
concept of linear normal modes to nonlinear systems. They have proved useful in a number 
of applications going from model reduction [1] and the study of localization phenomena [2] 
to system identification [3, 4], NNMs of undamped systems were first defined as a vibration 
of unison of the system [5]. The NNM definition was later generalized to encompass modal 
interactions [6]. Inspired by the center manifold technique, Shaw and Pierre extended the 
NNM concept to damped systems [7]. They defined NNMs as two-dimensional invariant 
manifolds in phase space which are tangent to the planes formed by the linear modes at the 
equilibrium point. Such manifolds are invariant under the flow (i.e. trajectories initialized 
in the manifold do not leave the manifold), a property which extends the invariance of linear 
modes to nonlinear systems. 
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Operator-theoretic methods have a long and rich tradition in dynamical systems. In 
particular, a (nonlinear) system can be studied through the spectral properties of the asso¬ 
ciated (linear) Koopman operator. Koopman operators were first defined for Hamiltonian 
systems in [8]. In the seminal work of Mezic [9] the spectral properties of the operator were 
connected to geometrical features of ergodic (i.e. conservative) systems. More recently, this 
spectral approach was extended to the case of dissipative systems in the context of model re¬ 
duction [10] and stability analysis [11]. The Koopman operator framework is also conducive 
to data analysis (see e.g. Dynamic Mode Decomposition [12]) and was successfully used to 
capture coherent structures (called Koopman modes) or detect instabilities in various fields 
(e.g. fluid flows [13], power grids [14], energy-efficient buildings [15]). 

Our aim is to explore the connections that exist between NNMs and the Koopman 
operator framework. In this direction, a fundamental result reported by Mezic in [9] is 
that normal modes of linear oscillations have their natural analogs, termed the Koopman 
modes, in the context of nonlinear dynamics. Specifically, it was shown that, for a hyperbolic 
fixed point, the dynamics can be linearized in the entire basin of attraction [16] and that 
the coordinates for the linearization can be obtained through the eigenfunctions of the 
Koopman operator [17]. This framework therefore seems very appropriate to extend the 
concept of linear normal modes to nonlinear systems. Indeed, we will demonstrate that 
there is a direct connection between NNMs and the spectral expansion of the Koopman 
operator. More precisely, the invariant manifolds are related to the zero level sets of some 
eigenfunctions of the operator, in direct analogy with the linear case. 

Using this new characterization, we will show how issues related to the parametrization of 
the invariant manifold can be overcome. In the original work of Shaw and Pierre, a single pair 
of variables (either in physical or modal coordinates) was used to obtain a parametrization of 
the manifold, which can possibly hold only locally. For instance, the parametrization may fail 
when the invariant manifold presents a complex geometry with foldings, as shown in [18, 19]. 
In contrast, interpreting the concept of NNMs in the Koopman operator framework allows 
us to obtain a parametrization that is valid globally. 

The paper is organized as follows. Sections 2 and 3 review the theory of NNMs and of the 
Koopman operator, respectively. In Section 4 the connection between the two frameworks is 
investigated, first in the linear case and then in the nonlinear case. We also derive a partial 
differential equation describing the geometry of the invariant manifold and an approximation 
of the manifold using Taylor series. The theoretical developments are illustrated with a two- 
degree-of-freedom mechanical system with cubic nonlinearity in Section 5. The conclusions 
of the present study are summarized in Section 6. 

2. Introduction to NNMs 

Consider the first-order equations of motion of an N-degree-of-freedom nonlinear me¬ 
chanical system 

Xi = Vi, yi = f(x, y) i = l,...,N (1) 

where aq represents a generalized coordinate (displacement or rotation), y t is the corre¬ 
sponding velocity, x = [aq, ...,aqv] and y = [yi, ...,yiy]. From this point on, we assume that 
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system (1) is damped, (x, y) = (0,0) is an equilibrium point, and the system linearized 
around (0,0) has n/2 pairs of complex-conjugate eigenvalues with non-zero real parts (hy¬ 
perbolic equilibrium point). 

In the 1990s, Shaw and Pierre defined a NNM as a two-dimensional manifold in phase 
space [7] that is invariant under the flow, i.e. trajectories with initial conditions in the 
manifold remain in it for all time. As such, the invariance property of linear normal modes 
was extended to nonlinear systems: a motion that starts on a mode stays on this mode, 
while the other modes remain quiescent. Inspired by the center manifold technique [20], 
Shaw and Pierre mathematically described the two-dimensional invariant manifold using a 
pair of master coordinates, the other state-space coordinates being functionally related to the 
master pair. Choosing arbitrarily (xj,i/j) as master coordinates, the remaining coordinates 
follow the constraint equations 

x i = Xj(xi,yi), yj = Y j (x i ,y i ), j = 1,... ,N ; j ± i. (2) 

In several other works (see, for instance, [21]), a linear change of coordinates was used to 
express Eq. (2) in modal space and describe the invariant manifold as a function of a linear 
modal coordinate and its corresponding modal velocity. 

A set of equations for the functions Xj and Y) can be derived by substituting the time 
derivative of Eq. (2) into the equations of motion (1). This leads to a set of 2N — 2 coupled 
partial differential equations (PDEs) 


w v dX j (x i ,y i ) dX j (x i ,y i ) 

_ <))', (./•/• //,) , ■ dYj ( x h yi ) 

Jj o IJi I Jil 

dxi dy.i 


with fj = fj(xi,X,yi,Y), X = {X j (x i ,y i )} 3 Y = {Y j (x i ,y i )}, and j = 1, ..., Y; j ± i. The 
PDEs (3) do not depend on time and can be solved for X and Y. Finally, after a substitution 
of the solutions in the ordinary differential equations governing the master coordinates, the 
dynamics on the NNM is reduced to a single-degree-of-freedom (SDOF) nonlinear oscillator 
dynamics 

x % = Vi, Vi = fi(xi,X,yi, Y). (4) 

The PDEs (3) admit several solutions that correspond to the extension of the underlying 
linear normal modes. The constraint relations (X, Y) describe the geometry of a single NNM 
in phase space. For a linear system, the mode corresponds to a plane whose dependence 
in Xi and y t is function of the inertia, stiffness and damping properties of the system. For 
a nonlinear system, the mode is no longer flat but curved due to nonlinear distortions. At 
the equilibrium point, the NNM is tangent to the plane formed by the linear mode of the 
underlying linear system. Figure 1 illustrates the in-phase NNM of a two-degree-of-freedom 
system with cubic nonlinearity. The surface shown in the figure corresponds to the position 
of the second DOF x 2 as a function of the master coordinates (xi,yi). The solutions of (3) 
can be approximated in the form of a polynomial series expansion in Xi and y t [7], whose 
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coefficients are the solutions of a set of algebraic equations. More recently, several numerical 
methods were proposed in the literature and applied to various conservative [21, 22] and 
nonconservative [19] systems, see [23] for a review. 



Figure 1: In-phase NNM of a two-degree-of-freedom system with cubic nonlinearity. The 
manifold was computed using the Taylor series expansion presented in Section 4.4 . 

Either derived using state space or modal coordinates, Eq. (2-4) are based on the fun¬ 
damental underlying assumption that the invariant manifold can be described as a function 
of the chosen pair of variables. However, the invariant manifold may present a complex 
geometry with several foldings that appear for a number of reasons including internal res¬ 
onances, localization, and multiple fixed points [23, 24, 25]. Such a manifold cannot gen¬ 
erally be described using a simple parametrization that is linear in the state space coor¬ 
dinates. For instance, the invariant manifold shown in Figure 1 exhibits foldings, which 
prevent a parametrization in terms of the master coordinates (xi,yi). In this case, a proper 
parametrization cannot be obtained as a solution of the PDEs (3). 

In Ref. [24], the parametrization issue was addressed using a larger set of master coordi¬ 
nates. When considering n m pairs, the slave coordinates can be described as 

•fij Aj (x m , y ni) j Uj y m ), (5) 

where (x m ,y m ) are the vectors of master coordinates. After solving the corresponding 
manifold-governing PDEs, the dynamics is reduced to m coupled (nonlinear) oscillators. 
Although this method is elegant and proved effective on a nonlinear beam [24, 25], it does 
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not completely solve the intrinsic parametrization issue and may therefore fail in some 
regions of the phase space. 

In this paper the parametrization issue is addressed using the relation that exists between 
NNMs and Koopman operator. Through a nonlinear change of coordinates, the dynamics 
on the NNM will be expressed as a singlc-degree-of-freedom linear oscillator. This novel 
approach bears strong resemblance with the normal form method [ 1 , 26], whose main idea 
is to simplify the vector held using a series of coordinate changes. Unlike the normal form 
method which has only local validity, the Koopman operator framework provides a change 
of coordinates that is valid in the entire basin of attraction of a hyperbolic fixed point [16]. 
Therefore, it is possible to think of methods based on Koopman operator as a generalization 
of methods based on normal form, with the notable difference that the system of coordinates 
obtained in this case is valid in the whole basin of attraction of a fixed point. By doing so 
the parametrization issue is overcome, and coupled PDEs describing the NNM in the entire 
basin of attraction of the equilibrium point can be obtained. 

3. Koopman Operator 

3.1. Motivation and definition 



Figure 2: 2DOF system with cubic nonlinearity. 

The system depicted in Figure 2 is considered in order to illustrate the Koopman operator 
definition and properties. The system dynamics can be expressed through a system of 
ordinary differential equations (ODEs) 

X = f(x), X = (xi, x 2 , x 3 , x a ) t ; ( 6 ) 

where x A and x 2 are, respectively, the positions of the first and second mass, and x 3 and 
x 4 the corresponding velocities. The general solution, or flow, of ( 6 ) is a function, <^>(x 0 ,f), 
that assigns a point in state space to an initial condition x 0 and a time t. 

It is often interesting to analyze the dynamics of a system in terms of specific quantities 
such as, for instance, the relative displacement between the first and second mass (?i(x) = 
x 2 — Xi, the potential energy ,r/ 2 (x) = \{k a (x\ + xf) + kf,(x i — x 2 ) 2 ) + \k n ix\ or the kinetic 
energy g 3 (x) = |(m 1 X 3 + m 2 x 2 ). All these quantities are functions from the state space to 
M (or more generally C) and are called observables. 
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Since state-space variables change with time, so do these observables. Given initial 
conditions x 0 , the time evolution of an observable is obtained through its composition with 
the flow. For instance, the potential energy varies according to V(x 0 ,t) = ^(^(xo, t)), while 
the kinetic energy is given by Jl(x 0 ,f) = ^(^(x 0 ,f)). At every time instant t, K(x. 0 ,t) 

and K(x 0 ,f) can still be considered as observables (i.e. functions from the state space to 

C) such that, generally speaking, the time evolution of an observable results in a family of 
observables given by its composition with the flow. 

The Koopman operator theory concentrates on how arbitrary observables change with 
time, instead of how trajectories are organized in state space. Given a dynamical system 

x = f(x), x G M n , (7) 

and a vector space of observables J 7 , the Koopman operator (or more rigorously, semigroup 
of operators) U l is defined as the composition of the flow with g(x) G J-\ 

U t g(x) = g(<l>(x.,t)). (8) 

Equation ( 8 ) considers the case of a scalar observable g. The definition can be extended 
to a vector-valued function g(x) by applying U t to each component g;(x) of the vector, 
provided that each vector component g?(x) is in T. 

The fundamental property of Koopman operator is linearity. Recalling the example in 
Figure 2, the evolution of the total energy of the system, < 74 (x) = < 72 (x) + < 73 (x), is given by 
the sum of the evolutions of g 2 and r/ 3 , i.e. E(x. 0 ,t) = K (x 0 ,f) + K(x 0 ,f). More generally, 
applying U l to the linear combination of two observables g Q (x) and gp (x) that belong to J- 
yields 


9 1 (x) = a g a (x) + b gp (x) (9) 

^Vr( x ) = 07 (^(M)) = ag a (<f>{x,t)) + b gp((j>(x,t)) = 

= atfgaix) + bU t gp(x). 

The linearity property of the Koopman operator is general and does not depend on the 
dynamical system. This property is very attractive, especially in the context of vibration 
analysis, where eigenvalues and eigenvectors form the basis of vibration modes. The Koop¬ 
man operator allows such tools to be used with nonlinear systems, but the price to pay is 
that infinite dimensional vector spaces have to be considered. 

3.2. Spectral properties 

From this point on the space of observables T is fixed. Although the only requirement 
is invariance with respect to the action of U t (given g G JF U t g G T for every t ), different 
choices of T result in different spectral properties. In this work, T is taken as the space 
of C 1 functions that are analytic in a neighborhood of the origin. Moreover, only systems 
with a stable hyperbolic fixed point (assumed to be the origin) are considered, and our 
investigations are restricted to the basin of attraction of that point. 
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Based on the above assumptions, this section presents the Koopman operator spectral 
properties that are essential to this paper. A more comprehensive introduction can be found 
in Ref. [17]. 

3.2.1. Eigenfunctions 

Since U 1 is linear, it is possible to consider its eigenfunctions: those observables s(x) G J- 
whose evolution in time is given by 

U t s(x ) = s(0(t,x)) = exp(At)s(x) Vt, (11) 

in which A G C is the corresponding eigenvalue. 

As an example consider a simple one-dimensional system 


x = —rjx, rj > 0, (12) 

with corresponding flow ( j>(x,t ) = xexp (—rjt). Applying the Koopman operator to the 
observable si(x) = x yields 

U t si(x) = si((f>(x,t)) = a; exp (—rjt) = exp(—rjt)si(x), (13) 


thus, the identity function si(x) = x is an eigenfunction with eigenvalue — rj. 

Now consider a general linear system defined by a matrix A. The eigenvalues of A are 
also part of the spectrum of the Koopman operator. Moreover, this property extends to 
nonlinear systems: if A is the Jacobian matrix evaluated at the origin of the vector held 
defined in (7), that is 


A 


df 



(14) 


then its eigenvalues A& are part of the spectrum of the Koopman operator [11], In fact, if A 
admits n linearly independent eigenvectors, its eigenvalues are associated with n independent 
eigenfunctions Sfc(x) that are defined in the entire basin of attraction of the equilibrium. 

However, the spectrum of U 1 is much richer than the spectrum of A: it always contains 
an infinity of eigenvalues. For instance, in the case of system (12), the observable s 2 (x) = x 2 
is also an eigenfunction: 


U t S 2 (x) = (4>(x,t)) 2 = (xexp(— gt)) 2 = x 2 exp(— 2ijt) = exp(—2 rjt)s 2 (x). (15) 

A similar reasoning can be applied to each integer power n of x, so that —nr) and s n (x) = 
s”(a;) = x n are eigenvalues and eigenfunctions of system (12), respectively. 

It turns out that this construction can be applied to general systems. Let Si(x) and 
S2 (x) be two eigenfunctions: 


fAsi(x) = exp(AR)si(x) 
TAs 2 (x) = exp(A 2 t)s 2 (x). 


(16) 
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Applying the Koopman operator to their product Si j2 (x) = Si(x)s 2 (x) yields 


U*s i >2 (x) = si >2 (^(x,t)) = Si(0(x,f))s 2 (0(x,f)) = 

= exp(Aif)si(x) exp(A 2 f)s 2 (x) = exp((Ai + A 2 )t)si(x)s 2 (x). 

Thus, the sum of two eigenvalues is also an eigenvalue, and equivalently the product of two 
eigenfunctions is an eigenfunction. Using (17) it is possible to construct an entire system of 
eigenfunctions. Considering the basis eigenfunctions Sk associated with the eigenvalues A/ c 
of the Jacobian matrix, we can obtain an infinity of eigenfuctions and eigenvalues: 

Sfci,...,fc n (x) = ^ :i ( x ) • • • s£ n (x), A kl ,...,kn = Mi + • • • + KK- (18) 

In the following a particular property of the eigenfunctions will be central to establish 
a connection to NNMs defined in Section 2. Let s(x) be an eigenfunction and A the cor¬ 
responding eigenvalue; then, its zero level set {x G M n j s(x) = 0 } is invariant for the 
dynamics (7) [11]. Indeed, given an initial condition xq satisfying 


s(x 0 ) = 0 , 


(19) 


the corresponding trajectory satisfies 

s(0(f,x o )) = exp(Af)s(x 0 ) = 0. (20) 

If the eigenvalue associated with s(x) is real, then so is s(x), and (19) corresponds to one 
equality. If A has a nonzero imaginary part, s(x) cannot be purely real, so that (19) is 
a complex equation yielding two equalities Re(s(x)) = 0 and Im(s(x)) = 0. In the latter 
case, the complex conjugate A is also an eigenvalue with eigenfunction s*,(x), so that (19) 
corresponds to the zero level set of two (complex conjugate) eigenfunctions. 

3.2.2. Spectral expansion 

One of the main tools in the analysis of linear systems is the expansion of the dynamics 
in terms of linear modes. For a general linear system 

X = Ax, X e M n , (21) 

with A diagonalizable, the flow can be written as 

n 

0(x,f) = 53s fc (x)v fc exp(Afct), ( 22 ) 

k=1 

where is a (possibly complex) eigenvector of A, A*, is the corresponding eigenvalue, and 
s*;(x) is the projection of x onto the eigenspace generated by v*,. Note the use of the 
same notation for the projections of x and the Koopman operator eigenfunctions; these two 
concepts will be linked in Section 4.1. 



The Koopman operator theory allows the decomposition (22) to be generalized to nonlin¬ 
ear systems [17]. Let g be an observable that belongs to the span of the eigenfunctions (18). 
The observable can be expressed as a (possibly infinite) sum of eigenfunctions 

3( x ) = Z C'fci,...,fc„4' 1 ( x )---^ n (x). (23) 

k\-\-...-\-k n ^>0 

Applying U 1 to (23), the time evolution of the observable is given by 

fA^x) = JZ C'fci,...,fc„«i 1 ( x ) • • • 4"( x ) exp((fciAi H-f k n X n )t). (24) 

(fci,...,fc„)eN n 

ki~\-...-\-kn >0 

The coefficients Ck lr ..,k„ are the so-called Koopman modes. Each of them is the projection of 
the observable g onto the eigenspace generated by the corresponding eigenfunction ... s^ n . 

If a vector-valued observable g : M n —> C d is considered, the Koopman modes in (23-24) 
become vector-valued modes Ck lr ..,k n and the scalar eigenfunctions remain unchanged. 

As an example, consider for system (12) the observable g±(x) = x + x 2 . In the previous 
section it was shown that the functions s k (x) = x k are eigenfunctions. Since gi(x) is a finite 
sum of eigenfunctions, we directly find the corresponding Koopman modes: C\ — 1, C 2 — 1, 

Ck = 0, k > 3. A less trivial example of spectral decomposition is given by the function 
gi(x) = exp (a;) — 1. In this case, the Taylor series 

OO-I OO 1 

g 2 (x) = exp (a;) -1 = 5] j-x k = J2 ( 25 ) 

k =1 K ' k= 1 K ~ 

shows that the Koopman modes are Ck — In both cases (24) can be used to obtain the 
time evolution of the observables. 

It is clear that (24) bears strong resemblance with the linear mode expansion (22). To 
obtain the same decomposition of the flow as in (22), it is necessary to find an observable 
whose evolution corresponds to the flow. This is achieved using the identity function id(x) = 
x. In a neighborhood of the origin, id(x) can be decomposed as sum of eigenfunctions [10] 

id(x) = x= 5] v fcl) ..., fcn sj 1 ( x )...s^(x). (26) 

(ki,...,k n )eN n 

k\-\- ...-\-kn >0 

Using (26) together with (24), a decomposition of the flow in terms of eigenfunctions is 
obtained 

lAid(x) = id (&(x,t)) =<£(£, x) = 5] v fcli ... ifcfl sj 1 (x)... s£ n (x) exp((&iAi + ... + k n X n )t). 

k±-\-.. .-\-k n ^>0 

(27) 

The decomposition (27) generalizes (22) to nonlinear systems. This infinite series can be 
truncated to a finite number of terms to approximate the system trajectories. The Koopman 
modes Vki,...,k„ associated with the identity function can be computed by solving a sequence of 
linear systems; the details of the calculation are reported in appendix A. For linear systems, 
Eq. (27) coincides with the linear normal mode decomposition, as shown in Section 4.1. 
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4. Koopman Operator and Vibration Modes 

J h l. The linear case 

Consider the linear system 

x = Ax, x G M n , (28) 

with A diagonalizable, and the corresponding expansion in linear modes given by 

n 

0(x,t) = Y Sfc(x)v fc exp(A fc t). (29) 

k =1 

The terms s&(x) are the projections onto the corresponding eigenspace; they can be written 
in the form 

Sfc(x) = w*x, (30) 

where * denotes the conjugate transpose and w/,. are the eigenvectors of A*, normalized so 
that 

* r 1 ■ h / q 1 \ 

w ‘ v ' 1 = { 0, k + h ■ (31) 

Applying the Koopman operator to Sfc(x), using the previous property, yields 

n 

f/ 4 s fc (x) = Sfc(<A(x, t)) = w* k Y Sj(x)vj exp(A^) = s fc (x) exp(A fc t). (32) 

i =1 

Thus, the projections Sfc(x) are eigenfunctions of U* (which explains the notation used 
for them). It follows that the expansions of the flow through linear modes and Koopman 
operator eigenfunctions are identical in the linear case. Also, the general decomposition (27) 
appears as the natural extension of linear normal mode decomposition to nonlinear systems. 

The relation between vibration modes and spectral expansion is well-known in the linear 
case: given a pair of complex conjugate eigenvalues Ai and A 2 , the motion on the linear 
mode is obtained when (29) reduces to 

si(x)vi exp(Ai t) + s 2 (x)v 2 exp(A 2 t), (33) 

i.e. when the following conditions on the eigenfunctions are verified: 

Sfc(x) = 0, k — 3,... ,n. (34) 

Thus, the plane defining the linear mode is determined by Eq. (34). Given a pair of complex 
conjugate eigenvalues, the corresponding linear mode is the intersection of the zero level set 
of all the other eigenfunctions. 

In the next two sections (34) is generalized to nonlinear systems, and it is shown that 
this generalization is equivalent to define NNMs in terms of invariant manifolds. 
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4-2. Nonlinear coordinate transformation 

The properties of Koopman operator in the general case of a nonlinear system with hy¬ 
perbolic fixed point can be studied through linearization techniques [16, 11]. The results 
of [16] are particularly useful, as they allow a local smooth change of coordinates (or diffeo- 
morphism) h(x) to be extended to the entire basin of attraction of the fixed point. Assuming 
that the eigenvalues are nonresonant (i.e. relations of the type rqA?: = 0 with n l e Id 
and J2i n i > 0 are not possible), there exists in a neighborhood of the origin an analytic 
diffeomorphism of the form [27] 

<9Z 

z = x + Z(x), Z(0) = 0, — 

such that 

z = Az. (36) 

Note that normal form methods aim at approximating the inverse transformation x = 
z + X(z) through its Taylor series. 

According to the results and methods given in [16], there exists a linearizing global C 1 
diffeomorphism h(x), which is analytic when restricted to a neighborhood of the origin and 
characterizes the eigenfunctions of the Koopman operator. In particular, the eigenfunction 
corresponding to the eigenvalue A/ c of A is given by 



»*(x) = wjh(x), 


(37) 


where w*. is a left eigenvector of A (here and in the following A is assumed diagonalizable). 
Moreover, it follows from (35) that 


<9h 

<9x 


x=0 


(38) 


so that 


ds k 

dx. 


= w 


ki 


x=0 


(39) 


i.e. the gradient of the eigenfunction at the origin is a left eigenvector of A. Note that 
the eigenfunctions s*,(x) = w^h(x) are related to the coordinates z. More precisely, they 
correspond to coordinates £ given by 


£ = Wz, 

( w* \ 

W= : , 

V < 

i.e. they are linearizing coordinates, in which the dynamics is diagonal. 


( 40 ) 
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4-3. The nonlinear case 

This section generalizes the connection between the eigenfunctions of Koopman operator 
and linear normal modes of vibration to NNMs as defined by Shaw and Pierre [7]. Recall 
from Section 2 that NNM manifolds are characterized by the following properties: 

1 . the manifold is invariant; 

2 . it passes through a stable fixed point (the origin); 

3. the plane tangent to it at the origin is a linear mode of the linearized dynamics. 

Following the discussion of Section 4.1, consider the surface defined by the intersection 
of n — 2 zero level sets of eigenfunctions 

Sfc(x) = 0, VA; ± k u k 2 , X kl = X k2 . (41) 

Relying on the properties of the previous sections, we can show that the manifold defined 
by (41) verifies the properties 1-3: 

1 . As shown in Section 3.2.1 the level sets s k (x) = 0 are all invariant, therefore the 
intersection of any number of them is invariant as well. 

2. From (35) and (37) it follows that 

Sfc(O) = w*h(0) = 0, (42) 

thus any zero level set contains the origin, and so does their intersection. 

3. The tangent space at the origin of (41) is defined by the linear equations 

f) Qi 

^(0)x = 0, (43) 

which, using (39), reduces to 

w£x = 0, V7c 7 ^ ki, k 2 , (44) 

and corresponds to the linear mode associated with and \ k2 as shown in Section 4.1. 

Therefore, the generalization of linear modes to nonlinear systems in the Koopman operator 
framework corresponds to the definition of NNMs proposed by Shaw and Pierre. 

4-4- NNM approximation using Taylor series 

Thanks to the connection established in the previous section, it is possible to use methods 
from the Koopman operator framework to compute NNMs. In this section a Taylor series 
approximation based on the spectral expansion (27) is proposed. The underlying idea is to 
use the values of the eigenfunctions as coordinates. This has the advantage that the dynamics 
in these coordinates is linear as illustrated in Section 4.2; furthermore, equations (41) have 
a particularly simple form in this coordinate system. 
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Let £& = Sfc(x) so that (26) becomes 

id(x)=x= £ v„. ■■■&’, (45) 

k\-\-.. .-\-kn^>0 

and (41) 

£* = Sfc(x) = 0 , Vfc ^ fci, k 2r X kl = X k2 . (46) 

Suppose ki — 1 and k- 2 = 2. Equation (45) expresses the transformation between the 
coordinate system defined by the Koopman operator eigenfunctions £ and the original co¬ 
ordinates x, while (46) defines the NNM in the £ coordinates. Combining the two, we can 
express the points on the NNM as a function of the variables £i and £ 2 : 

x= E v ^i fel £ 2 2 :=*( 6 , 6 ), (47) 

(fci,fc 2 )eN 2 

ki+k2>0 

where the notation v kl ,k 2 — v fci,fc 2 ,o,...,o is used. Equation (47) assigns to any pair of values of 
the eigenfunctions Si(x) and s 2 ( x ) a point on the manifold, and each point on the manifold 
corresponds to a value £i = £2 £ C. In other words, the manifold is parametrized using all 
the admissible values for the eigenfunctions. Moreover, it follows from (27) that the time 
evolution of a trajectory on the manifold is given by 

X W = E exp((fciAi + k 2 \ 2 )t) = tf(£i(t),£ 2 t(t)). (48) 

(fci,fc 2 )eN 2 

ki+k 2 >0 

Equation (47) can be truncated to a fixed order to approximate the invariant manifold 
in a neighborhood of the origin, and correspondingly (48) to approximate the NNM motion. 
The vectors v fc| k2 are the Koopman modes of the identity function and can be computed as 
discussed in Appendix A. 

To conclude this section some intrinsic limitations of approximating NNMs with the 
above Taylor series are discussed. First of all, the eigenfunctions might be analytic only 
in a neighborhood of the origin, so that the amplitude to which NNMs can be computed 
using (47) might also be limited (see also [11] for a similar issue in approximating eigen¬ 
functions with Taylor series). Even in cases in which the series converges in the whole basin 
of attraction, the computation of the Koopman modes is more demanding as the order 
increases, and becomes unfeasible for very high orders. 

Due to this the method proposed in this section is not meant to be a ready-to-use 
algorithm for the computation of NNMs of arbitrary systems, but rather a simple method 
to compute NNMs for limited amplitudes in small-dimensional systems. Its main purpose is 
to show that the theory of Koopman operator can lead to computational methods that use 
the coordinates given by the eigenfunctions. More generic algorithms could be derived using 
the results of the following section, in which a system of PDEs resembling (3) is derived. 
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4-5. Analogy with Shaw and Pierre's manifold governing PDE 

In Section 2 a system of PDEs (3) is derived starting from Shaw and Pierre definition of 
NNMs. A similar system can be derived for ^( 8 ) 8 ) defined in (47). To this end, note that 
the dynamics of 8 is by definition 

8 = Afc 8 - (49) 

It follows that the time derivative of \&(£i(t), 8 ( 8 ) is given by 


±(t) 


dti 


+ 


<98 


A 2 8(8- 


(50) 


Substituting (50) into the equations of motion (7), and evaluating at t — 0 yields 


m 


Wi 


A18 + 


<98 


A 2 8- 


(51) 


Equation (51) is a set of n PDEs very similar in principle to the one derived in Section 2. 
The approaches followed here and in Section 2 can be both interpreted as a way to solve 


Sfc(x) = 0, i = 3,...,N. (52) 

When using state-space variables as in (3), Eq. (52) is inverted, but the inversion might be 
valid only locally. In contrast, a different set of coordinates is used in (51) to parametrize the 
solution, which avoids issues arising when NNMs have a complex geometry (e.g. folding). 
The two systems of PDEs (3) and (52) have a similar structure, but (3) corresponds to 
71—2 equations while (52) consists of n equalities. Moreover, the second system admits 
more solutions, due to the fact that the eigenfunctions are defined modulo a scaling factor. 
Indeed, if VP( 8 , 8 ) is a solution then ^( 8 , 8 ) — ^( c 8 >c 8 ) is a solution for any complex 
c 7 ^ 0 . 

To obtain \l/( 8 , 8 ) real-valued 8 = 8 m| ist be assumed in (51). Equivalently, a system 
of equations can be derived using as variables u = Re( 8 ) and v = Im( 8 ), whose dynamics 

is 


u = <7i u — u±v, (53) 

v = u>iu + (Tiv, (54) 

with Ai = o\ + iuj i- This leads to 

w>= v ^(": r)(“)’ <55) 

with v ) = + iv, u — iv). Clearly, the above procedure can be followed for any pair 

of complex conjugate eigenfunctions, yielding a system of PDEs (51) or (55) for each NNM. 
These systems differ only for the value of A (or a and u). 
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4-6. Uniqueness of NNMs 

The section is concluded by a remark about conditions 1-3 of Section 4.3 defining the 
NNMs. Generally these conditions do not determine a unique manifold in the dissipative 
case. For the sake of illustration, consider a two-dimensional linear system 


x\ = ofiafi, 
x 2 = cr 2X2, 

with (j 2 < oq < 0. The phase portrait, illustrated in Figure 3, shows that almost all 
trajectories approach the linear eigenspace x 2 = 0, leading to a continuum of tangent mani¬ 
folds. Analogous phenomena are possible in higher dimensions. For completeness, Appendix 
B presents a family of invariant manifolds tangent to a linear normal mode, for a four¬ 
dimensional linear system. However, assuming nonresonant eigenvalues, the coefficients of 
a polynomial approximation are uniquely determined at any order, so that only one of the 
manifolds is analytic. This observation applies to Shaw and Pierre’s approach as well (note 
that the existence of locally analytic manifolds is guaranteed by the existence of a local 
analytic change of coordinates, see Section 4.2). The issue of non-uniqueness is discussed in 
detail in [28]. 

We note that also in the context of the Koopman operator there is the possibility to 
obtain multiple manifolds for only one eigenvalue. Specifically, our construction of the 
eigenfunctions relics on the analyticity of the transformation h of Section 4.2, which is 
uniquely defined by its Taylor series. However, if this smoothness condition is dropped it 
is possible to find different transformations verifying (35) and (36), obtaining in this way 
multiple eigenfunctions for the same eigenvalue. 


(56) 



xi 


Figure 3: Phase portrait of a typical two-dimensional linear system 
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5. A 2 DOF example 


The Taylor series expansion derived in Section 4.4 and based on the Koopman operator 
framework is used to compute the NNMs of the system shown in Figure 2. The equations 
of motion are 

miXx + k a x i + cxi + k nt x\ + k b (x i - x 2 ) + c(iq - x 2 ) = 0 , 
m 2 x 2 + k a x 2 + cx 2 + k b (x 2 - a+) + c(x 2 - iq) = 0, 

where mi = m 2 — rn — 1 Kg, c = 0.05 Ns/m, k a = 1 N/m and k n i = 0.5 N/m 3 . With 
yi = xi, y 2 = x 2 , ? 7 i = (k a + k b )/m , r] 2 = k b /m , f3 = c/m, a = k n i/m , (57) is recast in 
first-order form 


/ X 1 \ 


f 

0 

0 

1 

° \ 

l ) 


( 

0 

\ 

X 2 



0 

0 

0 

1 

x 2 

+ 


0 


Vi 



~Vi 

V‘2 

-2/3 

0 

Vi 


— axl 


V 2/2 / 


{ 

V2 

-Vl 

0 

- 2/3 ; 

V 2/2 ) 


\ 

0 

/ 


(58) 


The Jacobian matrix 


A = 


/ 0 
0 


0 
0 

Vi r l‘2 

V b2 -Vl 


1 

0 

-2/3 

P 


0 \ 

1 

/3 

-2/3 / 


(59) 


has two pairs of complex conjugate eigenvalues Ai = A 2 , A 3 = A 4 , which verify a 3 : 1 
resonance condition Ai = 3A 3 for k b = 4. Exact linear resonance is eliminated by varying k b 
but folding is likely to occur in the neighborhood of those parameter values. Table 1 shows 
how the eigenvalues vary for different values of k b . 


h 

Ai 

fn, 1 

Cl 

a 2 

to 

C 2 

4.7 

-0.075+13.224 

3.224 

0.0233 

-0.025+10.999 

1 

0.025 

4.3 

-0.075+13.079 

3.098 

0.0242 

-0.025+10.999 

1 

0.025 

4.1 

-0.075+13.032 

3.033 

0.0247 

-0.025+10.999 

1 

0.025 


Table 1: Eigenvalues, natural frequencies and damping ratios for different values of k b . 

To determine the coefficients describing the NNM in Eq. (47), the trajectory on the 
manifold given by (48) is injected into the equations of motion. Evaluating at time t — 0 
yields 


E V fcl,fc2^1 1 ^2 2 (^l^l + ^2 A 2 ) — 

(fci,fc 2 )eN 2 

ki+k2>0 


16 



( 


0 

0 


(60) 


\ 


A E v fel)fc2 e^ 

(fcl,fc 2 )6N 2 

ki+k 2 >0 


V 


tt 

> 

\ 

a 

E V fcl 


V 

1 (fci,fc 2 )eN 2 
\ fci+fc 2 >0 

) 

1 ) 


o 


/ 


) 


where (v) i stands for the ith component of the vector v. Matching the terms corresponding 
to the first-order powers results in two systems 


(Ail — A)v 10 — 0, (A 2 I — A)v 0 i — 0, 


(61) 


so that the first-order modes v 10 and v 0 i are the eigenvectors of A corresponding to Ai and 
A 2 . If there is no resonance between the eigenvalues of A, the successive systems admit a 
unique solution, and are obtained by matching like-power terms (see the appendix for the 
complete derivation). Doing so for the second-order terms yields linear equations of the form 

((AqAi + &2A2)I — A)vfc 1 fc 2 = 0, k\ + k 2 = 2, (62) 


where the right hand side is zero due to the absence of second-order terms in the vector 
held. We can proceed iteratively, solving for any mode a linear system of the form 

(rl - A)v fclfc2 = b, r = £qAi + k 2 \ 2 (63) 

where b is obtained as a sum of lower-order modes (see appendix A for the details). A 
change of coordinates using the Jordan basis is used to improve numerical stability. 

The above procedure is used to compute the high-frequency (i.e. the out-of-phase) mode 
of the 2 DOF system for kb = 4.3. The Taylor series used to generate the figures of this 
section were truncated at the 50 th order. The maximum amplitude at which the manifolds 
are represented was chosen so that trajectories obtained with (48) agree with those computed 
through numerical integration. To have a quantitative measure of the difference between 
two trajectories the normalized mean square error (NMSE) was used: 

1 DO K 

NMSE= — y ||x(fc)-x(fc)|| 2 , (64) 

U X k =0 

where x is the time series obtained with (48), x is the reference series obtained through 
numerical integration, and crj: i s ^ ie variance of the latter. The NMSEs obtained to validate 
the manifolds were always lower than 1%, a value usually taken as threshold for a good 
match between two time series. 

The result is compared with the finite element (FE) method developed in [19]. This 
numerical method is based on an explicit parametrization of the invariant surface using a 
pair of master coordinates and it is known to be accurate provided that a good mesh is 
employed. For the out-of-phase mode, the selected master coordinates correspond to the 
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displacement and velocity of the first mass. As can be seen in Figure 4 the results from the 
two methods agree very well, but the FE method allows higher levels of amplitude to be 
reached, while the Taylor series validity is more limited. 

The lower-frequency (in-phase) NNM was computed for the different values of kb reported 
in Table 1. The results are shown in Figure 5 and compared with the NNMs obtained using 
the FE algorithm. As kb is decreased from 4.7 (Figures 5a and 5b) to 4.3 (Figures 5c 
and 5d) and 4.1 (Figures 5e and 5f), the manifold geometry changes and starts to fold, 
which prevents the use of master coordinates for the parametrization. This is also confirmed 
by Figure 6, where the manifold obtained for kb = 4.1 is shown with a viewpoint orthogonal 
to the X 2 I /2 plane. The FE algorithm reaches the folding but cannot go beyond, since it 
cannot converge in region where the master coordinates parametrization is not valid. In 
contrast the folding does not pose any problem with the parametrization obtained through 
the Koopman eigenfunctions. 

In Figure 7 one of the trajectories considered for the validation is compared with the 
corresponding numerical integration, in the case kb = 4.1. As expected according to the 
values of the NMSE, the two trajectories cannot be distinguished by visual inspection. 

Equation (48) also highlights a structure of the evolution that is reminiscent of that 
obtained using Fourier series with periodic trajectories. The main difference is that the 
number of terms in (48) is not the same at all orders, but increases with the order. Figure 8 
shows the terms corresponding to the Erst three non-zero orders, for kb = 4.7 and kb = 
4.1, and in Figure 9 their sum is compared with the trajectories obtained using numerical 
integration. When only terms up to order five are used the agreement between the two 
trajectories is not as good as in Figure 7, nonetheless, it is remarkable that the main aspects 
of the dynamics are retained by this approximation, at least in the case considered herein. It 
is also interesting to note that in the case kb = 4.7 the trajectory presents a lower harmonic 
content, and using the same order a better approximation is obtained. 
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Figure 5: In-phase NNM for different values of k b ; (a)-(b) k b = 4.7, (c)-(d) k b = 4.5, (e)-(f) 
k b = 4.1;^ FE method, i=i Taylor series 19 
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(a) (b) 


Figure 6: Different views of the invariant manifold shown in Figure 5 (A;& = 4.1 ); (a) X\ 
(Figure 5e) top view, (b) X 2 (Figure 5f) bottom view. Note that the manifold undergoes 
folding, which prevents a parametrization with state-space variables. 
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(a) 



(c) 


Figure 7: Comparison of time series;— 
initial condition: ^i(O) = —0.23, £ 2 ( 0 ) = 
^ 2 , (c) 2 / 1 , (d) y 2 



t 


(b) 



(d) 


Numerical integration,— Taylor approximation; 
-0.76, 2/i(0) = -0.4, 2 , 2 (0) = -0.86; (a) m, (b) 
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X\ 



Figure 8: Decomposition of trajectories in different terms;—1 st order term,—3 rd order 

term,. 5 th order term; (a) kb = 4.7 (xi(0) = —0.57, x 2 (0) = —0.55, 2/1 (0) = —0.44, 

2/2(0) = -0.66); (b) k b = 4.1 (xi(0) = -0.23, x 2 (0) = -0.76, yi (0) = -0.4, y 2 ( 0) = -0.86) 



Figure 9: Comparison of 5 th order approximation and numerical integration;— Numerical 
integration,— 5 th order approximation; (a) kb = 4.7 (xi(0) = —0.57, x 2 (0) = —0.55, 2 / 1 ( 0 ) = 
-0.44, 2/ 2 (0) = -0.66), (b) k b = 4.1 (:n(0) = -0.23, x 2 (0) = -0.76, 2 / 1 ( 0 ) = -0.4, 2 / 2 ( 0 ) = 
- 0 . 86 ) 
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6. Conclusion 


This research investigates the intimate connection that exists between NNMs of clamped 
systems and the Koopman operator associated with their dynamics. In particular, we have 
shown that there is a correspondence between the invariant manifolds defined by Shaw and 
Pierre and the zero level sets of particular eigenfunctions of the Koopman operator, and 
that this correspondence is a natural generalization of the linear case. 

Through the Kooman operator eigenfunctions a new characterization of NNMs is possi¬ 
ble, and a novel approach to their parametrization is proposed. The main advantage of this 
new approach is its global validity, obtained as a consequence of the linearizing properties 
provided by the eigenfunctions. This parametrization overcomes a known limitation in the 
use of master coordinates. 

However, the Taylor series expansion used herein is local in nature and restricted to 
low-dimensional systems. In order to target arbitrarily high levels of amplitude, as well 
as high-dimensional systems, more efficient methods should be developed for solving the 
system of PDEs (51) or (55). As mentioned in Section 4.5 there are differences between 
these systems and the one obtained in the work of Shaw and Pierre. Nonetheless, all these 
systems share a common structure, which might allow the extension of numerical methods 
developed using master coordinates to this new framework. In addition, the Koopman 
operator approach is associated with recent data analysis techniques that are well-suited to 
high-dimensional systems [29]. This could potentially open the door to develop data-driven 
methods for computing NNMs of large systems. 

The preliminary results presented here motivate a deeper investigation of NNMs for 
damped systems. For example it is well known that in conservative systems NNMs can 
bifurcate, and exceed the number of degrees of freedom of the system, but how these phe¬ 
nomena survive the presence of damping is still an open question. 
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A. Computation of the Koopman Modes of the Identity 

In this appendix we derive a method for computing the Koopman modes of the 

identity (see (26)). In particular, we show that these modes are the solutions of linear 
systems. 

We fix some notation. For a = (an, ..., a n ) £ N n let 

n 

a! = Y[ aj\, (65) 

3 =1 
22 



( 66 ) 


a 


d Q {, = 


E a ji 

3 =1 

ai“if i 


dxi 1 ... dx% n ’ 


= n x ? 

3 =1 


Using this notation, we can write 

Sfc lr ..,fc„(x) =s fc (x), fc = (/ci,...,/c n ), 

with 


s(x) = 


In addition, (26) can be rewritten as 


/ si(x) > 

V s ^(x) ) 


X = E E v fc s fc (x), 


and 


h>1 fceN™ 
|fc|=/l 


</>(*, x ) = f ( x ) = E E v fc s fc (x)/cA, 

/i>lfceN n 


where k\ = YJj=\ kj\j and A j are the eigenvalues of the Jacobian matrix 


A = 


df 

<9x 


x=0 


Similarly, the Taylor series of f(x) can be written as 

fiW = £ £ 


h>1 fceN n 


and injecting (71) in (74) leads to 


$(*)=£ £ 


hi>l k i€N" 
\ki\=h! 


h\ 


3 =1 


\ 




E E ( v fc 2 )t5 fc2 ( x ) 


/i2>i feeN" 
\ 1^21=^2 


J 


By comparison of (72) and (75), for any k G N n we obtain 


T, 3 s 


Mik \=eee n n k. 

i=l 

23 


h-l 


(67) 

( 68 ) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 

(75) 

(76) 



These equations are linear in the terms with maximum value of \k\, and can be solved if all 
the Vfc/ with \k'\ < \k\ are known. 

The matrix defining the linear system related to has the form kXI — A. The first-order 
modes are the eigenvectors of A and it is necessary to choose a scaling factor, which cor¬ 
respondingly fixes the value of the eigenfunctions (defined up to a multiplicative constant). 
If there are no resonances between the eigenvalues, the linear systems related to modes of 
higher order are uniquely determined, since they are defined by a nonsingular matrix. We 
can also use the known form of the matrices defining the linear systems to simplify the com¬ 
putation. In the case of a diagonalizablc Jacobian matrix, the systems are diagonal when 
expressed in coordinates defined by a Jordan basis. 

To obtain a manifold using (47) not all the modes are necessary, but just a subset of them: 
the Vfc, k = (&!,..., k n ), for which only two indices kj 1 and kj 2 are nonzero. Since to obtain 
the linear system corresponding to a given mode v^, only the modes v^./, k! = (k[,..., k\ J, 
with ki > k{ are needed, the computation can be restricted to the modes that appear in (47). 

As concluding remarks we note that approximating a solution of (51) with a Taylor series 
yields the same linear systems. Also, the procedure presented here corresponds in essence to a 
normal form approximation. In that context, starting from a dynamical system with diagonal 
linear part, a transformation can be obtained so that the resulting equations of motion 
are the simplest possible. The only algebraic condition needed to obtain a transformation 
that linearizes the dynamic is the nonresonant eigenvalues. The corresponding linearizing 
coordinates are the eigenfunctions of the Koopman operator. 


B. A family of invariant manifolds for a linear system 

Consider a four-dimensional linear system 

x = Ax, (77) 

with 

—G\ 0J\ 0 0 

—uq —0i 0 0 

0 0 — 0" 2 0J 2 

0 0 — u 2 —cr 2 

Note that any four-dimensional linear system with two couples of complex conjugate eigen¬ 
values can be transformed into this form. It will be assumed that 0 \ > 0 and a 2 > 0 
(stability). 

For any value of ci and c 2 consider the manifold x(u, v ) = f (u, v; c i, c 2 ), 




fl(u, 

v) 

= u 


f2(u, 

v) 

= V 





a 2 

h (u, 

V) 

= (V + u 2 ; 

| 2o-f 




cr 2 

h(u, 

V) 

= [u 2 + V 2 ] 

1 2cr i 



/log (u 2 + v 2 ) w 2 \ 

V 2^ ) 

/log (u 2 + v 2 ) u 2 \ 

V 2oq ) 

24 


Ci — sin 
Ci + cos 


/log [u 2 + u 2 )o; 2 \ 

V 2oq ) 

/log (u 2 + u 2 )o; 2 \ 

V 2oq J 



( 79 ) 



The derivatives of f are 


fuA = 1 
fu.2 = 0 


(72 _^ 

u(u 2 + V 2 ) 2rT i (cos 


fu, 3 


log(« 2 +t> 2 )u >2 

2(71 


(cicr 2 - c 2 co 2 ) - sin 


logfu 2 +v 2 )u>2 

2(71 


(c 2 a 2 + ciu 2 ) 


fuA = 


<?1 

u (u 2 + V 2 )^~ l fcosf^_±Ti^) (c 2 (7 2 + CiU 2 ) + sill ) (ci<r 2 - c 2 u; 2 ) 


o-i 


i = 0 

fv, 2 = 1 


fv.'.i ~ 


fvA = 


v (tt 2 + v 2 ) 2ct i 1 ^cos (^ log ( u 2 + 0 ( Ci(72 : _ c 2 u; 2 ) _ S i n j d° 6 (^ 2 +» )^ ( C2(72 + Cl o; 2 ) 


o-i 


V (u 2 + v 2 ) 2ct i 1 ^cos ^ lQg ^ 2 ^ — ) (c 2 cr 2 + Ci^ 2 ) + sin ( Cl<j2 _ P^ 2 ) 


Cl 


(80) 


These are continuous for a 2 > cy, and more regular if the distance between the two real 
parts of the eigenvalues is increased further. 

Assuming sufficient regularity, f defines a manifold, with tangent plane at a point (u, v ) 
spanned by the vectors 


df 


df 


f — _ f — 

L U 1 L V r\ ' 

ou ov 


(81) 


Since it holds 


Af (u,v) = ( — <7iU + UJ\v)f u + (— UJ\U — <7iv)i v , 


(82) 


all the manifolds are invariant ((82) is the analogous of the equation derived by Shaw and 
Pierre equation for the system (77)). Moreover, if a 2 > a i 


f«(0,0) 


( l \ 


f °\ 

0 

0 

, f„(0,0) = 

1 

0 

V 0 ) 


1 0 / 


(83) 


so that at the origin the tangent plane is 0:3 = 0 , 2:4 = 0 , which corresponds to the linear 
normal mode, but only for ci = c 2 = 0 the whole manifold is the linear normal mode. 

Note that apart from the points verifying x lj0 = a: 2 ,o = 0 (i.e. points on the fast linear 
mode), for any other point (x li0 , x 2)0 , x 3 ) 0 , 2 : 4 , 0 ), there exists a manifold which passes through 
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it, this is obtained using the parameters 



As a final remark we note that even if each couple of coefficients c±, C 2 defines a different 
manifold, the derivatives of all the functions defined by (79) have the same values at the 
origin (up to the order at which they are defined): they all share the same polynomial 
approximation. Interestingly the same situation is encountered in the context of center 
manifolds. 
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